Нам понадобятся следующие пакеты.
c-CS-s: 9 mice c-CS-m: 10 mice c-SC-s: 9 mice c-SC-m: 10 mice
t-CS-s: 7 mice t-CS-m: 9 mice t-SC-s: 9 mice t-SC-m: 9 mice
d <- d[,-1]
aggr(d, sortVars=F, combined=T, bars=F, numbers=T, prop=F, sortCombs=T) #see NAs
Мы видим 3 наблюдения с NA в больше половине случаев, их я удалю. Всего 552 наблюдения с полной информацией.
d <- d[- (which(rowSums(is.na(d))>10)),] ##Exclude obs with many NAs
print("Количество оставшихся NA")
## [1] "Количество оставшихся NA"
sort(colSums(is.na(d)),decreasing=T)[1:10]
## BCL2_N H3MeK4_N BAD_N EGR1_N H3AcK18_N pCFOS_N ELK_N
## 285 270 213 210 180 75 15
## Bcatenin_N MEK_N DYRK1A_N
## 15 4 0
Посмотрим на pointrange.
d_d <- d
d_d$class <- as.factor(d_d$class)
ggplot(d_d, aes(x = d_d$class, y = d_d$BDNF_N, colour = class)) +
labs(y = 'BDNF_N', x = 'Class') +
theme_light() +
stat_summary(geom = 'pointrange', fun.data = mean_cl_normal)
По графику видно, что группы из середины отличаются. Сделаем дисперсионный анализ
# функция Anova из пакета car
mod_1 <- lm(BDNF_N ~ class, data = d_d)
m_anova <- Anova(mod_1)
# встроенная функция
m1 <- aov(d$BDNF_N ~ d$class, data = d_d)
summary(m1)
## Df Sum Sq Mean Sq F value Pr(>F)
## d$class 7 0.2878 0.04112 18.82 <2e-16 ***
## Residuals 1069 2.3362 0.00219
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# просто сравнила функции - показывают одно и тоже :)
По данным дисперсионного анализа, мы видим, что уровень BDNF_N достоверно зависит от класса животного.
# Данные для графиков остатков
mod_diag <- fortify(mod_1)
# График расстояний Кука
ggplot(mod_diag, aes(x = 1:nrow(mod_diag), y = .cooksd)) +
geom_bar(stat = 'identity')
Остатки распределнены нормально.
ggplot(mod_diag, aes(x = mod_diag$class, y = .stdresid)) + geom_boxplot()
Остатки по классам не различаются.
qqPlot(mod_1, id = FALSE)
Распределение очень близко к нормальному.
Сделаем пост-хок тест, чтобы узнать, какие группы различаются. Я взяла тест Тюкей, поскольку он относительно не строгий и не слабый.
posthoc <- glht(mod_1, linfct = mcp(class = "Tukey"))
summary(posthoc)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = BDNF_N ~ class, data = d_d)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## c-CS-s - c-CS-m == 0 0.0030979 0.0055459 0.559 0.99930
## c-SC-m - c-CS-m == 0 -0.0482717 0.0053980 -8.942 < 0.001 ***
## c-SC-s - c-CS-m == 0 -0.0258249 0.0055459 -4.657 < 0.001 ***
## t-CS-m - c-CS-m == 0 -0.0264852 0.0055459 -4.776 < 0.001 ***
## t-CS-s - c-CS-m == 0 -0.0337570 0.0059483 -5.675 < 0.001 ***
## t-SC-m - c-CS-m == 0 -0.0181541 0.0055459 -3.273 0.02392 *
## t-SC-s - c-CS-m == 0 -0.0136310 0.0055790 -2.443 0.22105
## c-SC-m - c-CS-s == 0 -0.0513696 0.0055459 -9.263 < 0.001 ***
## c-SC-s - c-CS-s == 0 -0.0289228 0.0056900 -5.083 < 0.001 ***
## t-CS-m - c-CS-s == 0 -0.0295831 0.0056900 -5.199 < 0.001 ***
## t-CS-s - c-CS-s == 0 -0.0368549 0.0060829 -6.059 < 0.001 ***
## t-SC-m - c-CS-s == 0 -0.0212520 0.0056900 -3.735 0.00474 **
## t-SC-s - c-CS-s == 0 -0.0167289 0.0057223 -2.923 0.06894 .
## c-SC-s - c-SC-m == 0 0.0224468 0.0055459 4.047 0.00155 **
## t-CS-m - c-SC-m == 0 0.0217865 0.0055459 3.928 0.00226 **
## t-CS-s - c-SC-m == 0 0.0145147 0.0059483 2.440 0.22285
## t-SC-m - c-SC-m == 0 0.0301176 0.0055459 5.431 < 0.001 ***
## t-SC-s - c-SC-m == 0 0.0346406 0.0055790 6.209 < 0.001 ***
## t-CS-m - c-SC-s == 0 -0.0006603 0.0056900 -0.116 1.00000
## t-CS-s - c-SC-s == 0 -0.0079321 0.0060829 -1.304 0.89717
## t-SC-m - c-SC-s == 0 0.0076708 0.0056900 1.348 0.87968
## t-SC-s - c-SC-s == 0 0.0121939 0.0057223 2.131 0.39494
## t-CS-s - t-CS-m == 0 -0.0072718 0.0060829 -1.195 0.93310
## t-SC-m - t-CS-m == 0 0.0083311 0.0056900 1.464 0.82595
## t-SC-s - t-CS-m == 0 0.0128542 0.0057223 2.246 0.32461
## t-SC-m - t-CS-s == 0 0.0156029 0.0060829 2.565 0.16969
## t-SC-s - t-CS-s == 0 0.0201260 0.0061130 3.292 0.02262 *
## t-SC-s - t-SC-m == 0 0.0045231 0.0057223 0.790 0.99360
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
Построим график, чтобы визуализировать результат (одинаковые буквы означают неразличающиеся группы)
MyData <- data.frame(class = factor(levels(d_d$class), levels = levels(d_d$class)))
MyData <- data.frame(MyData, predict(mod_1, newdata = MyData, interval = 'confidence'))
ggplot(data = MyData, aes(x = class, y = fit)) +
geom_bar(stat = 'identity', aes(fill = class), width = 0.5) +
geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.1)+
labs(x = 'Class', y = 'Экспрессия белка BDNF_N') +
geom_text(aes(y = 1.6, label = c('A', 'A', 'B', 'C', 'D', 'E', 'A', 'A')), colour = 'black', size = 5, position = position_stack(vjust = 0.5))+
theme(legend.position = 'none')
Поскольку у нас около 75 белков, до построения модели, попробуем сократить количество анализируемых переменных. ### 1. Анализ на мультиколлинеарность. Попарно сравним коэффициент корреляции белков.
m <- length(colnames(d))
# Функция для расчета p-value для коэффициента корреляции (взято: http://www.sthda.com/english/wiki/visualize-correlation-matrix-using-correlogram)
# mat : is a matrix of data
# ... : further arguments to pass to the native R cor.test function
cor.mtest <- function(mat, ...) {
mat <- as.matrix(mat)
n <- ncol(mat)
p.mat<- matrix(NA, n, n)
diag(p.mat) <- 0
for (i in 1:(n - 1)) {
for (j in (i + 1):n) {
tmp <- cor.test(mat[, i], mat[, j], ...)
p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
}
}
colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
p.mat
}
# matrix of the p-value of the correlation
d <- as.data.frame(d)
p.mat <- cor.mtest(d[,-c((m-3):m)])
head(p.mat[, 1:5])
## DYRK1A_N ITSN1_N BDNF_N NR1_N NR2A_N
## DYRK1A_N 0.000000e+00 0.000000e+00 3.436576e-34 5.382704e-23 6.337133e-28
## ITSN1_N 0.000000e+00 0.000000e+00 1.798784e-57 7.246317e-48 9.222316e-49
## BDNF_N 3.436576e-34 1.798784e-57 0.000000e+00 7.780844e-247 3.484596e-195
## NR1_N 5.382704e-23 7.246317e-48 7.780844e-247 0.000000e+00 0.000000e+00
## NR2A_N 6.337133e-28 9.222316e-49 3.484596e-195 0.000000e+00 0.000000e+00
## pAKT_N 2.215088e-09 1.111704e-06 1.177053e-26 2.320727e-12 2.918397e-04
cor <- cor(d[,-c((m-3):m)],use = "pairwise.complete.obs")
corrplot(cor, type = 'upper', p.mat = p.mat,
sig.level = 0.01, insig = 'blank',
tl.col="black", tl.srt=45,
tl.cex = 0.4)
Данных очень много, поэтому крайне высока вероятность наличия мультиколлинеарности. Одним из признаков мультиколлениарности является высокая парная зависимость между малозначимыми переменными, поэтому найдем пары белков с наибольшей коллинеарностью и удалим одного из пары.
## [1] "Number: Correlations between proteins > 0.8"
## [1] 87
## [1] "Number: Correlations between proteins > 0.5"
## [1] 435
87 белков из датасета имеют корреляцию больше 0,8. Удалим каждого из пары таких белков.
Var_na <- which(colSums(is.na(d))>0)
d <- as.data.frame(d)
cor_na <- cor(d[,-c(Var_na,(m-3):m)],d[,Var_na],use = "pairwise.complete.obs")
d <- d[, which(! colnames(d) %in% names(which(colSums(cor_na > 0.8) > 0)))]
m <- length(colnames(d))
Vars_Mcol <- findCorrelation(cor(d[,-((m-3):m)],
use = "pairwise.complete.obs"),0.8)
d <- d[,- Vars_Mcol]
Осталось 54 белков.
m <- length(colnames(d))
d <- as.data.frame(d)
aov <- numeric()
for( i in (1: (m-4))){
aov[i] <- summary(aov(d[,i]~d$ERBB4_N))[[1]]$"Pr(>F)"[1]
}
print("number of protein variables which have a significant relationship to ERBB4_N")
## [1] "number of protein variables which have a significant relationship to ERBB4_N"
sum( aov < 0.05/length(aov))
## [1] 37
Vars <- which(aov < 0.05/length(aov))
d <- d[,c(Vars,(m-3):m)]
Итого осталось 37 белков, которые значимо коррелируют с ERBB4_N. Для построения обычной линейной модели предикторов слишком много, сделаем РСА.
which(colnames(d)=='ERBB4_N')
## [1] 23
m <- length(colnames(d))
# уберем колонки c не белками
prot_only <- as.matrix(d[,c(-1,-23, -(m-3):-m)])
# стандартизируем переменные
# Function to standardize numeric variables
z_numeric_vars<-function(data_frame) {
for (n in names(data_frame)) {
if (class(data_frame[[n]]) == "numeric" | class(data_frame[[n]]) == "integer"){
var = paste(n,"_z", sep="")
data_frame[[var]] <- scale(data_frame[[n]], center = TRUE, scale = TRUE)
data_frame[[n]] = NULL
}
}
data_frame
}
prot_st <- z_numeric_vars(prot_only)
# заменим NA на среднее
prot <- na.aggregate(prot_st)
Построим ординацию
Построим график собственных чисел
screeplot(ord, bstick = TRUE, type = 'lines')
## [1] 52.5942
Оставим первые 3 компоненты, они объясняют 53% изменчивости.
Построим факторные нагрузки
biplot(ord, scaling = 'species', correlation = TRUE,
main = 'PCA - species scaling', display = 'species')
Построим график для первых трех компонент.
df_scores <- data.frame(d, scores(ord, display = 'sites', choices = c(1, 2, 3), scaling = 'sites'))
plot_ly(df_scores, x = ~PC1, y = ~PC2, z = ~PC3, color = df_scores$class, size = 0.5)
prot_1 <- as.data.frame(prot)
ERBB4_N <- scale(d$ERBB4_N)
new_d <- as.data.frame(cbind(ERBB4_N, scores(ord, display = 'sites', choices = c(1, 2, 3), scaling = 'sites')))
lmod <- lm(V1 ~ ., data = new_d)
summary(lmod)
##
## Call:
## lm(formula = V1 ~ ., data = new_d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.54187 -0.36359 -0.02816 0.34886 2.60972
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.829e-16 1.941e-02 0.000 1
## PC1 3.365e+00 9.421e-02 35.716 < 2e-16 ***
## PC2 -7.180e-01 1.107e-01 -6.486 1.35e-10 ***
## PC3 2.133e+00 1.323e-01 16.125 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6371 on 1073 degrees of freedom
## Multiple R-squared: 0.5952, Adjusted R-squared: 0.5941
## F-statistic: 525.9 on 3 and 1073 DF, p-value: < 2.2e-16
Модель объясняет 59% изменчивости.
Проверка на мультиколлинеарность.
sqrt(vif(lmod))>2
## PC1 PC2 PC3
## FALSE FALSE FALSE
Мультиколлинеарности нет.
Однако, vif выдаёт такой странный результат. (если можно, поясни, пожалуйста, в каком месте я ошиблась или что с этим делать)
vif(lmod)
## PC1 PC2 PC3
## 1 1 1
Проверим остатки на нормальность
qqPlot(lmod, labels = row.names(states), simulate = TRUE, main = 'График Q-Q')
## sit142 sit250
## 142 250
Есть какие-то выбросы, но в целом похоже на нормальное распределение.
Проверим остатки на линейность.
crPlots((lmod))
Если нелинейные связи есть, то они не очень большие.
Построим визуализацию.
new_data <- data.frame(ERBB4_N = seq(min(new_d$V1), max(new_d$V1), length.out = 50),
PC1 = mean(new_d$PC1),
PC2 = mean(new_d$PC2),
PC3 = mean(new_d$PC3))
new_data$predicted <- predict(lmod, newdata = new_data) # предсказанные значения
new_data$SE <- predict(lmod, newdata = new_data, se.fit = TRUE)$se.fit # стандартные ошибки
new_data$upper <- new_data$predicted + 2 * new_data$SE # верхняя граница доверительной области
new_data$lower <- new_data$predicted - 2 * new_data$SE # нижняя граница доверительной области
Визуализируем влияние первой компоненты.
ggplot(data = new_data, aes(x = ERBB4_N, y = predicted)) +
geom_ribbon(data = new_data, aes(ymin = lower, ymax = upper), alpha = 0.3) +
geom_line(color = 'blue', size = 1) +
geom_point(data = new_d, aes(x = ERBB4_N, y = PC1)) +
labs(y = 'PC1', x = 'ERBB4_N')
Построить модель не получилось, облако точек явно меняется по другой закономерности.
m <- length(colnames(d))
# уберем колонки c не белками
prot_only <- as.matrix(d[,c(-1, -(m-3):-m)])
# стандартизируем переменные
# Function to standardize numeric variables
z_numeric_vars<-function(data_frame) {
for (n in names(data_frame)) {
if (class(data_frame[[n]]) == "numeric" | class(data_frame[[n]]) == "integer"){
var = paste(n,"_z", sep="")
data_frame[[var]] <- scale(data_frame[[n]], center = TRUE, scale = TRUE)
data_frame[[n]] = NULL
}
}
data_frame
}
prot_st <- z_numeric_vars(prot_only)
# заменим NA на среднее
prot <- na.aggregate(prot_st)
Построим ординацию
Построим график собственных чисел
screeplot(ord, bstick = TRUE, type = 'lines')
## [1] 52.89198
Оставим первые 3 компоненты, они объясняют 53% изменчивости.
Построим факторные нагрузки
biplot(ord, scaling = 'species', correlation = TRUE,
main = 'PCA - species scaling', display = 'species')
Построим график для первых трех компонент.
df_scores <- data.frame(d, scores(ord, display = 'sites', choices = c(1, 2, 3), scaling = 'sites'))
plot_ly(df_scores, x = ~PC1, y = ~PC2, z = ~PC3, color = df_scores$class, size = 0.5)